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Abstract 

In this paper, we present a study of an a posteriori estimator for the 
discretization error of a non-standard finite difference scheme applied to 
boundary value problems defined on an infinife inferval. In particular, we 
show how Richardson’s exfrapolafion can be used fo improve fhe numerical 
solution involving fhe order of accuracy and numerical solutions from fwo 
nested quasi-uniform grids. A benchmark problem is examined for which 
fhe exacf solution is known and we gel fhe following resull: if fhe round-off 
error is negligible and fhe grids are sufficienlly fine Ihen fhe Richai'dson’s 
error eslimale gives an upper bound of fhe global error. 
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1 Introduction 


The main aim of this paper is to show how Richardson’s extrapolation can be used 
to define an error estimator for a non-standard finite difference scheme applied to 
boundary value problems (BVPs) defined on the infinite interval. Without loss of 
generality, we consider the class of BVPs 


^ = f(x,u), xe[0,oo), 

g(u(0),u(°o)) = 0, 


( 1 . 1 ) 


where u(x) is a <7—dimensional vector with ^u{x) fox as components, 

f : [0, oo) X —)■ R^, and g : R^ x R^ —)■ R^. Here, and in the following, we use 
Lambert’s notation for the vector components [15, pp. 1-5]. Existence and unique¬ 
ness results, as well as results concerning the solution asymptotic behaviour, for 
classes of problems belonging to (1.1) have been reported in the literature, see 
for instance Granas et al. [12], Countryman and Kannan [7], and Agarwal et al. 
[2,3, 1]. 

Numerical methods for problems belonging to (1.1) can be classified accord¬ 
ing to the numerical treatment of the boundary conditions imposed at infinity. 
The oldest and simplest treatment is to replace the infinity with a suitable finite 
value, the so-called truncated boundary. However, being the simplest approach 
this has revealed within the decades some drawbacks that suggest not to apply 
it especially if we have to face a given problem without any clue on its solution 
behaviour. Several other treatments have been proposed in the literature to over¬ 
come the shortcomings of the truncated boundary approach. In this research area 
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they are worth of eonsideration: the formulation of so-ealled asymptotie bound¬ 
ary conditions by de Hoog and Weiss [9], Lentini and Keller [16] and Markowich 
[17, 18]; the reformulation of the given problem in a bounded domain as studied 
first by de Hoog and Weiss and developed more recently by Kitzhofer et al. [14]; 
the free boundary formulation proposed by Fazio [10] where the unknown free 
boundary can be identified with a truncated boundary; the treatment of the orig¬ 
inal domain via pseudo-spectral collocation methods, see the book by Boyd [5] 
or the review by Shen and Wang [21] for more details on this topic; and, finally, 
a non-standard finite difference scheme on a quasi-uniform grid defined on the 
original domain by Fazio and Jannelli [11]. 

When solving a mathematical problem by numerical methods one of the main 
concerns is related to the evaluation of the global error. For instance, Skeel [23] 
reported on thirteen strategies to approximate the numerical error. Here we are 
interested to show how within Richardson’s extrapolation theory we can derive an 
error estimate. For any component U of the numerical solution, the global error e 
can be defined by 

e = u-U, (1.2) 

where u is the exact analytical solution component. Usually, we have several 
different sources of errors: discretization, round-off, iteration and programming 
errors. Discretization errors are due to the replacement of a continuous problem 
with a discrete one and the related error decreases by reducing the discretiza¬ 
tion parameters, enlarging the value of N, the number of grid points in our case. 
Round-off errors are due to the utilization of floating-point arithmetic to imple¬ 
ment the algorithms available to solve the discrete problem. This kind of error 
usually decreases by using higher precision arithmetic, double or, when available, 
quadruple precision. Iteration errors are due to stopping an iteration algorithm 
that is converging but only as the number of iterations goes to infinity. Of course, 
we can reduce this kind of error by requiring more restrictive termination criteria 
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for our iterations, the iterations of Newton’s method in the present ease. Program¬ 
ming errors are beyond the scope of this work, but they can be eliminated or at 
least reduced by adopting the so-called structured programming. When the nu¬ 
merical error is caused prevalently by the discretization error and in the case of 
smooth enough solutions the discretization error can be decomposed into a sum 
of powers of the inverse of N 


/lyo 

u — Un + GqI—J -{-Cl 



+ C2 



(1.3) 


where Q, Ci, C 2 , ... are coefficients that depend on u and its derivatives, but are 
independent on N, and po, pi, P 2 , ■■ ■ are the true orders of the error. The value of 
each for k = 0, 1, 2, • • •, is usually a positive integer with Pq < pi < pi < ■■ ■ 
and all together constitute an arithmetic progression of ratio pi — po, see Joyce 
[13]. The value of po is called the asymptotic order or the order of accuracy of the 
method or of the numerical solution U. 


2 Numerical scheme 

In order to solve a problem in the class (1.1) on the original domain we discuss 
first quasi-uniform grids maps from a reference finite domain and introduce on the 
original domain a non-standard finite difference scheme that allows us to impose 
the given boundary conditions exactly. 

2.1 Quasi-uniform grids 

Let us consider the smooth strict monotone quasi-uniform maps x = x{^), the so- 
called grid generating functions, see Boyd [5, pp. 325-326] or Canuto et al. [6, p. 
96], 

x=-c-ln{l-^) , (2.1) 
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and 


x = c 


1 -^ ’ 


( 2 . 2 ) 


where ^ G [0,1], a: G [0, oo], and c > 0 is a eontrol parameter. So that, a family of 
uniform grids = n/A^ defined on interval [0,1] generates one parameter family 
of quasi-uniform grids Xn = x{^n) on the interval [0,°°]. The two maps (2.1) and 
(2.2) are referred as logarithmie and algebraie map, respectively. As far as the 
authors knowledge is concerned, van de Vooren and Dijkstra [24] were the first 
to use this kind of maps. We notice that more than half of the intervals are in the 
domain with length approximately equal to c and xyv_i = clnA for (2.1), while 
~ cN for (2.2). For both maps, the equivalent mesh in x is nonuniform 
with the most rapid variation occurring with c <^x. The logarithmic map (2.1) 
gives slightly better resolution near x = 0 than the algebraic map (2.2), while the 
algebraic map gives much better resolution than the logarithmic map as x —)■ oo. In 
fact, it is easily verified that 


c-ln(l-<^) 


for all ^, but ^ =0, see figure 1 below. 

The problem under consideration can be discretized by introducing a uniform 
grid of A -f 1 nodes in [0,1] with ^q = 0 and = ^n + h with h=\ /N, so 
that Xn is a quasi-uniform grid in [0, oo]. The last interval in (2.1) and (2.2), namely 
[xAf_i,Viv], is infinite but the point X[^_i /2 is finite, because the non integer nodes 
are defined by 



with n G {0,1,..., A — 1} and 0 < a < 1. These maps allow us to describe the 
infinite domain by a finite number of intervals. The last node of such grid is placed 
at infinity so right boundary conditions are taken into account correctly. Figure 1 
shows the two quasi-uniform grids x = x„, n = 0,1,..., A defined by (2.1) and by 
(2.2) with c = 10 and A equal to, from top to bottom, 10, 20 and 40, respectively. 
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Figure 1: Quasi-uniform grids: top frame for (2.1) and bottom frame for (2.2). 
We notice that, in both cases, the last mesh-point is xn = °°- 


In order to derive the finite difference formulae, for the sake of simplicity, we 
consider a generic scalar variable u(x). We can approximate the values of this 
scalar variable at mid-points of the grid by 


^u+3/4 “^fi+l/2 , ^n+1/2 ~^«+l/4 

•^n4-3/4 '*^«-|-l/4 '*^11+3/4 ^n+l/A 


(2.3) 


As far as the first derivative is concerned we can apply the following approxima¬ 
tion 

Un+l U,i 

_ 

n-|-l/2 2 (x„_j_ 3^4 


du 

dx 
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(2.4) 






















































These formulae use the value um = u{°°), but not xm = For a system of differ¬ 
ential equations, (2.3) and (2.4) can be applied component-wise. 

2.2 A non-standard finite difference scheme 

A non-standard finite difference scheme on a quasi-uniform grid for the class of 
BVPs (1.1) can be defined by using the approximations given by (2.3) and (2.4) 
above. A finite difference scheme for (1.1) can be written as follows: 

— Un — «n+l/2t(-*^n-|-l/2;^;j+l/2Un+l +C„+l/2Un) = 0 , 

for n = 0,l,...,A-l (2.5) 

g(Uo,Uyv)=0, 


where 


^n+l/2 

^«+l/2 

1 /2 


2 (^«+3/4 “^n-l-1/4) ; 
•^«+l/2 “•^n+1/4 
•^n+3/4 “•^n+1/4 
•^«+3/4 —^n+l/2 

5 

''•n+3/4 


( 2 . 6 ) 


for n = 0,1,...,A—1. The finite difference formulation (2.5) has order of accu¬ 
racy 0{N~^). It is evident that (2.5) is a nonlinear system of d (A-|-1) equations in 
the (A^ + 1) unknowns U = (Uo,Ui,... ,UAf)^. For the solution of (2.5) we can 
apply the classical Newton’s method along with the simple termination criterion 


1 

d{N+l) 


d N 

££|A^f/„|<TOL, 


£= 1 n=0 


(2.7) 


where for n = 0,1,..., A and £ — 1,2,..., J, is the difference between two 
successive iterate components and TOL is a fixed tolerance. 
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3 Richardson’s extrapolation and error estimate 


The utilization of a quasi-uniform grid allows us to improve our numerieal results. 
The algorithm is based on Riehardson’s extrapolation, introdueed by Riehardson 
in [19, 20], and it is the same for many finite differenee methods: for numerieal 
differentiation or integration, solving systems of ordinary or partial differential 
equations, see, for instanee, [22]. To apply Riehardson’s extrapolation, we earry 
on several ealeulations on embedded uniform or quasi-uniform grids with total 
number of nodes Ng for g = 0, 1 ,..., G: e.g., for the numerieal results reported in 
the next seetion we have used 5,10,20,40, 80,160, 320, 640,1280,2560, or 5120 
grid-points. We ean identify these grids with the index g = 0, the eoarsest one, 
1, 2, and so on towards the finest grid denoted by g = G. Between two adjaeent 
grids, all nodes of largest steps are identieal to even nodes of denser grid due to the 
uniformity. To find a more aeeurate approximation we ean apply k Riehardson’s 
extrapolations on the used grids 


Ug+i,k+i 


Ug+i,k + 


^g+l,k Ug,k 

2Pk - 1 


(3.1) 


where g e {0,1,2,.. .,G — 1}, G {0,1,2,... ,G — 1}, 2 = Ng+i/Ng appearing in 
the denominator is the grid refinement ratio, and pk is the true order of the dis- 
eretization error. We notiee that to obtain eaeh value of f/^+i requires having 
eomputed two solution U in two embedded grids, namely g -f 1 and g at the ex¬ 
trapolation level k. For any g, the level k = 0 represents the numerieal solution 
of U without any extrapolation, whieh is obtained as deseribed in subseetion 2.2. 
In this ease, Richardson extrapolation uses two solutions on embedded refined 
grids to define a more accurate solution that is reliable only when the grids are 
sufficiently fine. The case k = I is the classical single Richardson’s extrapolation, 
which is usually used to estimate the discretization error or to improve the solution 
accuracy. If we have computed the numerical solution on G -|- 1 nested grids then 
we can apply equation (3.1) G times performing G Richardson’s extrapolations. 
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The theoretical orders pk of accuracy of the numerical solution U with k ex¬ 
trapolations verify the relation 


Pk = Po + k{Pi -Po) . 


(3.2) 


where this equation is valid for G {0,1,2,..., G — 1}. In any case, the values of 
Pi^ can be obtained a priori by using appropriate Taylor series or a posteriori by 


log( \Ug,k-^\)-'^ogi\Ug+i,k-u\) 

log(2) 


(3.3) 


where u is again the exact solution (or, if the exact solution is unknown, a reference 
solution computed with a suitable large value of N) evaluated at the same grid- 
points of the numerical solution. 

To show how Richardson’s extrapolation can be also used to get an error esti¬ 
mate for the computed numerical solution we use two numerical solutions and 
U 2 N computed by doubling the number of grid-points. Taking into account equa¬ 
tion (3.1) we can conclude that the error estimate by Richardson’s extrapolation is 
given by 


_ UlN — Un 
~ 2P0 - 1 


(3.4) 


where po is the true order of the discretization error. Hence, E is an estimation 
of the truncation error found without knowledge of the exact solution. We notice 
that E is the error estimate for the more accurate numerical solution U 2 N but only 
on the grid points of Ga?. 


4 Numerical results: a BVP in colloids theory 

In this section, we consider a benchmark problem with known exact solution for 
our error estimator. It should be mentioned that all numerical results reported in 
this paper were performed on an ASUS personal computer with iV quad-core Intel 
processor and 16 GB of RAM memory running Windows 8.1 operating system. 
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The non-standard finite difference scheme described above has been implemented 
in FORTRAN. The numerical results reported in this section were computed by 
setting 

T0L=1E-12. (4.1) 


The benchmark problem, see Alexander and Johnson [4], arises within the 
theory of colloids and is given by 


cp-u 

— 2sinh(M) = 0 xG[0,oo], 

m(0) = uq , u{°°) = 0 , 


(4.2) 


where uq > 0. The exact solution of the BVP (4.2) 


u{x) = 2 In 


(e“o/2 + l) + 

(^^o/2 -|- 1) ^V^^ — (^^o/2 — t) 


(4.3) 


has been found by Countryman and Kannan [7, 8], and the missing initial condi¬ 
tion is given by 

^(0) = -2 v^cosh(Mo) - 1 . (4.4) 

We rewrite the governing differential equation as a first order system and 
indicate the exact solution with u = and the numerical solution with 

U = (^[/, ^t/)^. In order to fix a specific problem, as a first test case, we consider 
Mo = 1- As mentioned before we used 5,10,20,40, 80,160, 320, 640,1280,2560, 
or 5120 grid-points, so that G = 10, and we adopted a continuation approach for 
the choice of the first iterate. This means that the accepted solution for A = 5 is 
used as first iterate for A = 10, where the new grid values are approximated by 
linear interpolations, and so on. The first iterate for the grid with A = 5, where the 
field variable was taken constant and equal to one and its derivative was taken also 
constant and equal to minus one, is shown in the top frame of figure 2. The bottom 
frame of the same figure shows the accepted numerical solution. Our relaxation 
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a: 



Figure 2: Sample iterates for problem (4.2) with uq = 1. 

algorithm takes seven iterations to verify the termination eriterion (2.7) with TOL 
given by (4.1). Onee the eontinuation approach has been initialized, the iteration 
routine needs 3 or 4 iterations to get a numerical solution that verifies the stop¬ 
ping criterion. For the sake of completeness in figure 3 we display the numerical 
solution for 77 = 40 along with the exact solution. 

Figure 4 shows in a log by log scale the computed errors. We can compute the 
order of accuracy po, p\ and p 2 according to the formula (3.3). As it is easily seen 
from figure 4 we got po ~ 2, 4 and /?2 ~ 6 for both the field variable and 

its first derivative. As far as the a posteriori error estimator is concerned in figure 
5 we report the computation related to two sample cases: namely, the estimate 
obtained by using N = 20, 40 and N = 40, 80. We notice that the global error, for 
both the solution components, is of order 10“^ and it decreases as we refine the 


11 















1^ 




0.5 


s 

(N 

s' 


0 


;:d 

-0.5 


0 



I s m « 


• * 


A 


A 


-1 - 


-1.5^ 


o 

‘u 


-V 

V 

'u 

A 



2 3 4 5 


Figure 3: Final iterate and exact solution for problem (4.2) with uq = \ and N = 
40.' zoom of the transitory region. 




Figure 4: Graphical derivation of the values ofpQ, p\ and P 2 for the non standard 
finite difference scheme applied to (4.2) with uq = 1. 

grid. It is easily seen that the estimator defined by equation (3.4) provides upper 
bounds for the global error. 
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Figure 5: Global and a posteriori error estimates for the field variable and its first 
derivative for (4.2) with uq = 1. Top: N = 20, 40, and bottom: N = 40, 80. Here 
^ e and ^e are the global errors by equation (1.2), whereas and are the error 
estimates provided by equation (3.4). 
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A more challenging test case is given by setting uq = 1 . In figure 6 we display 
the numerical solution for A = 5120 along with the exact solution. In table 1 



Figure 6: Final iterate and exact solution for problem (4.2) with uq = 1 and N = 
5120.' zoom of the transitory region. 

we list the computed as well as the extrapolated values obtained for the missing 
initial condition. For the sake of brevity, in this table, we do not report the fewer 
accurate values obtained with the coarser grids. These results can be compared 
with the exact value, ^(0) ~ —46.789615734913319, obtained by equation (4.4). 

Figure 7 shows in a log by log scale the computed errors. From figure 7 it is 
clear that the computed orders, using equation (3.3), are slightly different from 
the theoretical ones, namely: po ~ 1-99, pi 3.96 and p 2 ~ 5.77. 

As far as the a posteriori error estimator is concerned in figure 8 we report the 
computation related to two sample cases: namely, the estimate obtained by using 
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g,o 




gA 


2f/, 


g.2 


160 -43.835177171609345 

320 -45.864298511341850 -46.540672291252690 

640 -46.537797149336093 -46.762296695334179 -46.777071655606278 
1280 -46.725033491934731 -46.787445606134277 -46.789122200187620 

2560 -46.773360098843838 -46.789468967813541 -46.789603858592159 

5120 -46.785544794016836 -46.789606359074504 -46.789615518491907 
Table 1: Richardson’s extrapolation for ^(0) = ^Uq. 




Figure 7: Problem (4.2) with uq = 1, global errors for: the field variable on the 
left and its first derivative on the right. 


N = 1280, 2560 and N = 2560, 5120. Once again the global error, for both the 
solution components, decreases as we refine the grid and the estimator defined by 
equation (3.4) provides upper bounds for the global error. 
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Figure 8: Zoom in the domain related to the initial transitory for global errors 
and a posteriori error estimates for the field variable and its first derivative for 
(4.2) with uq = 1. Top: N = 1280, 2560 and bottom: N = 2560, 5120. Here ^e 
and are the global errors by equation (1.2), whereas and are the error 
estimates provided by equation (3.4). 
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5 Conclusions 


In this paper, we have defined a posteriori estimator for the global error of a non¬ 
standard finite differenee seheme applied to boundary value problems defined on 
an infinite interval. A test problem was examined for whieh the exaet solution 
is known and we tested our error estimator for two sample eases: a simpler one 
with smooth solutions eomponent and a more ehallenging one presenting an ini¬ 
tial fast transitory for one of the solution eomponents. For this seeond test ease, 
we showed how Riehardson extrapolation ean be used to improve the numeri- 
eal solution using the order of aeeuraey and numerieal solutions from two nested 
quasi-uniform grids. Moreover, the reported numerieal results elearly show that 
our non-standard finite differenee seheme, implemented along with the error esti¬ 
mator defined in this work, ean be used to solve ehallenging problems arising in 
the applied seienees. 

In our previous paper [11] we derived instead of equation (2.3) the finite dif¬ 
ferenee formula 

Un+l/2 ~ ^- Un H---M„+l . (5.1) 

However, by setting n = N this formula (5.1), replaeing xjy = reduees to 
“ 7 V- 1/2 = “iv-i that does not involve the boundary value Uf^ and therefore the 
boundary eondition eannot be used. In that work, this was the reason that foreed us 
to modify this formula at n = A — I, see [ 11 ] for details. The mentioned drawbaek 
is eompletely overeome by the new formula (2.3). 
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